rm(list = ls(all = TRUE)) #clear workspace

library(foreign)
library(survey)
options(survey.lonely.psu = "adjust")

setwd("C:/Users/martin.romero/Dropbox/ENVIPE/2013/BASES")
#setwd("~/Dropbox/ENVIPE/2013/BASES")

################################################
#Pregunta al entrevistado
################################################

data <- read.dbf("tper_vic.DBF")

#Ponderador y diseno
data$FAC_ELE <- as.numeric(as.character(data$FAC_ELE))
design.per <- svydesign(id=~UPM,strata=~EST, weights=~FAC_ELE, data=data)

#Lo secuestraron para exigir dinero o bienes:
#1 Sí
#2	No
#9	No sabe / no responde
table(data$AP7_3_12) #80 casos
foo <- svymean(~AP7_3_12, design=design.per, na.rm=TRUE) #0.11%
print(foo)
confint(foo)
foo <- svytotal(~AP7_3_12, design=design.per, na.rm=TRUE) #88,526
print(foo)
confint(foo)
cv(foo)

#7.4 Â¿Me podrÃ�a decir cuÃ¡ntas veces sufriÃ³ (CÃDIGO DE LA INCIDENCIA)?
table(as.numeric(as.character(data$AP7_4_12)))
sum(as.numeric(as.character(data$AP7_4_12)), na.rm=TRUE)
foo <- svytotal(~as.numeric(as.character(AP7_4_12)),design=design.per, na.rm=TRUE) #89,086
print(foo)
confint(foo)
cv(foo)

################################################
#Modulo de victimizacion
################################################

modulo <- read.dbf("tmod_vic.dbf")
    
#6.2 ¿Me podría decir si sus secuestradores.    				
#1  exigieron rescate a familiares por su liberación?
#2	lo obligaron a retirar dinero de un cajero, entregar joyas, celular u otras cosas?
#9	No sabe / no responde
#b	blanco
table(modulo$BP6_2)
round(prop.table(table(modulo$BP6_2)),3)*100

#6.1 ¿Cuánto tiempo lo tuvieron secuestrado?  					
#1  Menos de 24 horas
#2	De 1 a 3 días
#3	De 4 a 10 días
#4	De 11 a 29 días
#5	De 1 a 3 meses
#6	Más de 3 meses
#9	No sabe / no responde
#b	blanco
table(modulo$BP6_1)
round(prop.table(table(modulo$BP6_1)),3)*100

      
################################################
#
#E S T I M A D O R   P O R   H O G A R
#
################################################      

#Ponderador y diseno
data$FAC_HOG <- as.numeric(as.character(data$FAC_HOG))
design <- svydesign(id=~UPM,strata=~EST, weights=~FAC_HOG, data=data)

########################################
# Son integrantes del hogar?
########################################

#6.10a Â¿CuÃ¡ntos integrantes del hogar sufrieron secuestro?
table(data$AP6_10_2, useNA="ifany")
sum(as.numeric(as.character(data$AP6_10_2)), na.rm=TRUE) # 278 victimas de secuestro originalmente reportados
# estimacion de victimas sin correccion: 111,005
svytotal(as.numeric(as.character(data$AP6_10_2)), design=design, na.rm=TRUE)

#6.11 De los integrantes que me acaba de mencionar, Â¿cuÃ¡les vivÃ�an y compartÃ�an los alimentos en este hogar?
table(data$AP6_11_1, useNA="ifany")
table(data$AP6_11_2, useNA="ifany")
table(data$AP6_11_3, useNA="ifany")
#Notese que 43+10+1=53 personas secuestradas no son del hogar (sin contar cuarta persona)

es.hogar <- data[,c("AP6_11_1","AP6_11_2","AP6_11_3")]
es.hogar <- apply(es.hogar, 2, as.numeric)
es.hogar[es.hogar %in% c(NA,2)] <- 0
sum(es.hogar) #verificar que suma 221

######################
# Numero de secuestros
######################
#6.12 Â¿CuÃ¡ntas veces sufriÃ³ el secuestro el integrante del hogar (NÃMERO)?
secuestros <- cbind(data[,c("AP6_12_1","AP6_12_2","AP6_12_3")])
secuestros <- apply(secuestros, 2, as.numeric)
secuestros[is.na(secuestros)] <- 0
sum(secuestros)

#Verificar que pertenecen al hogar
secuestros <- cbind(victima.1=secuestros[,1]*es.hogar[,1],
                    victima.2=secuestros[,2]*es.hogar[,2],
                    victima.3=secuestros[,3]*es.hogar[,3])
  
#Sumar la cuarta victima del hogar (asumimos un secuestro)
table(data$AP6_10_2)
victima.4 <- data$AP6_10_2=="04"
victima.4[is.na(victima.4)] <- 0
secuestros <- cbind(secuestros,victima.4)
secuestros[rowSums(es.hogar)==0] <- 0 # hogar con missreport

# Numero de eventos en la encuesta
sum(secuestros) #239 secuestros reportados
svytotal(rowSums(secuestros), design=design) #Estimacion: 105,682

####################
# Numero de victimas
####################
victimas <-  rowSums(secuestros>0)
sum(victimas) #223 victimas reportadas
svytotal(victimas, design=design) #Estimacion: 94,438

####################
# Numero de hogares
####################
hogares <- victimas>0
table(hogares) #197 hogares con al menos un secuestro
svytotal(hogares, design=design) #Estimacion:84,487

#Estimacion original, sin corregir:
#6.10 Durante el 2012 en este paÃ�s (MÃ©xico), Â¿algÃºn integrante de este hogar 
#sufriÃ³ un secuestro o secuestro exprÃ©s, para exigir dinero o bienes?
table(data$AP6_10_1)
svytotal(~AP6_10_1, design=design) #Estimacion sin corregir: 95,280

################################################
#Intervalos de confianza
################################################
ci.data <- data.frame(Cluster=data$UPM, 
                      Strata=data$EST, 
                      Weights=data$FAC_HOG, 
                      secuestros=rowSums(secuestros), victimas, hogares)
ci.design <- svydesign(id=~Cluster,strata=~Strata, weights=~Weights, data=ci.data)

# Al 90%
confint(svytotal(~hogares, design=ci.design), level=0.9)
confint(svytotal(~victimas, design=ci.design), level=0.9)
confint(svytotal(~secuestros, design=ci.design), level=0.9)

# Al 95%
confint(svytotal(~hogares, design=ci.design), level=0.95)
confint(svytotal(~victimas, design=ci.design), level=0.95)
confint(svytotal(~secuestros, design=ci.design), level=0.95)

# Coeficientes de variacion
cv(svytotal(~hogares, design=ci.design))
cv(svytotal(~victimas, design=ci.design))
cv(svytotal(~secuestros, design=ci.design))

################################################
#Informantes que se equivocaron al definir hogar
################################################

error <- data[,c("AP6_11_1","AP6_11_2","AP6_11_3")]
error <- apply(error, 2, as.numeric)
error[error %in% c(NA,1)] <- 0
error <- rowSums(error)
sum(error>0) #45 informantes con error



